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Abstract 

Genetic variation within plant species is determined by a number of factors such as reproductive mode, breeding system, 
life history traits and climatic events. In alpine regions, plants experience heterogenic abiotic conditions that influence the 
population's genetic structure. The aim of this study was to investigate the genetic structure and phylogeographic history 
of the subalpine shrub Rhododendron ferrugineum across the Pyrenees and the links between the populations in the 
Pyrenees, the Alps and Jura Mountains. We used 27 microsatellite markers to genotype 645 samples from 29 Pyrenean 
populations, three from the Alps and one from the Jura Mountains. These data were used to estimate population genetics 
statistics such as allelic richness, observed heterozygosity, expected heterozygosity, fixation index, inbreeding coefficient 
and number of migrants. Genetic diversity was found to be higher in the Alps than in the Pyrenees suggesting colonization 
waves from the Alps to the Pyrenees. Two separate genetic lineages were found in both the Alps and Pyrenees, with a 
substructure of five genetic clusters in the Pyrenees where a loss of genetic diversity was noted. The strong differentiation 
among clusters is maintained by low gene flow across populations. Moreover, some populations showed higher genetic 
diversity than others and presented rare alleles that may indicate the presence of alpine refugia. Two lineages of R. 
ferrugineum have colonized the Pyrenees from the Alps. Then, during glaciation events R. ferrugineum survived in the 
Pyrenees in different refugia such as lowland refugia at the eastern part of the chain and nunataks at high elevations leading 
to a clustered genetic pattern. 
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Introduction 

The maintenance of genetic diversity within populations is 
essential for the conservation of species since it can enable 
adaptation to future environmental changes [1], [2]. With high 
genetic diversity, it is more likely that some individuals will possess 
alleles that could allow the population to adapt to a new habitat 
[3]. High genetic diversity may be maintained via gene flow 
among individuals and among populations. In recent years, many 
studies have shown strong impacts of global change on species' 
genetic diversity [4]. Because climate warming is greater at high 
than at low elevations [5] alpine biodiversity hotspots that are 
home to many endemic and rare species may be particularly 
threatened by climate change [6]. Thus, to predict how these 
biodiversity hotspots will deal with climate changes, it is important 
to understand the genetic patterns of the species and their changes 
over time [7]. 

Genetic structure and genetic divergence within and among 
plant populations result from a number of contemporary and 
historical factors acting at various temporal and spatial scales. 
Alpine habitats are typically characterized by small spatial changes 
in environmental variables such as topography, temperature, snow 
cover, soil moisture or bedrock [8] favorable to phenotypic 
differentiation and local genetic variability over short geographical 



distances. In addition to abiotic factors, life history traits can also 
influence genetic features of populations. For instance, selfmg 
perennials and clonal species generally exhibit lower levels of 
genetic diversity and higher levels of population differentiation 
than outcrossing perennials [9] . In contrast, mixed mating systems 
maintain a high level of genetic diversity and may be advantageous 
for plants in variable alpine environments [6]. The combined 
effects of both life history traits and small-scale abiotic heteroge- 
neity may result in local genetic adaptations and high genetic 
differentiation of populations, notably where gene exchanges 
between populations are restricted as is often the case in mountain 
habitats [10]. Therefore, plant-community fragmentation pattern, 
steep climatic gradients leading to asynchronous flowering 
phenology or disjoint pollinator communities [7], and topograph- 
ical barriers may hamper gene flow (through seed and pollen 
movements) and may increase differentiation between populations. 

At large scales, species genetic patterns are often shaped by past 
climate-driven range dynamics [11]. In temperate regions, species 
ranges have changed dramatically over the Pleistocene with its 
successions of glaciations and warming periods. In mountains, 
periodic advances and retreats of glaciers have, on several 
occasions, forced species to descend to lowland and peripheral 
refugia and/or to survive in isolated ice-free areas above glaciers 
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(nunataks). Genetic isolation, bottlenecks and potential allopatric 
evolution of populations in refugia have been followed by repeated 
founding events during population expansions, which have led to 
loss of genetic diversity and heterozygosity [11]. Additionally, 
populations located at eroding margins of species distribution 
ranges following climate warming could have persisted in suitable 
isolated lowland habitats [12]. It has been suggested [13-14] that 
such small populations isolated for so long could display low 
within-population genetic diversity (as a result of genetic drift) but 
have disproportionally high levels of between-population differen- 
tiation and regional genetic diversity and distinctiveness [15]. 

The genetic pattern and phylogeography of many species [16], 
[17] have been investigated in various mountainous regions, but to 
date, only a very limited number of studies have focused on the 
Pyrenees [18-20]. The Pyrenees are a mountain range stretching 
across the isthmus that lies between the Iberian Peninsula and the 
rest of the European continent. Acting as a biogeographical 
barrier during postglacial expansions, they are a limit for 
numerous Iberian and northern species. However, lower lands at 
the eastern and western extremities of the mountain range have 
been corridors for migration of several plant and animal species 
[21,22]. Due to oceanic influences to the west and Mediterranean 
influences to the east of the chain, the climate varies greatly along 
a longitudinal gradient in the Pyrenees (Meteo France data). This 
bioclimatic pattern potentially allowed some species to survive 
glaciations in various refugia along the chain. Studies have 
suggested that species have survived glaciations [23-25] in two 
main refugia - one in south central and another in the eastern 
Pyrenees - thanks to mild Mediterranean and oceanic climatic 
influences. However both geological and biological data [26] 
demonstrate that during glacial periods the ice sheet did not 
completely cover the highlands and rarely descended below 
1000 m a.s.l except for big glaciers pushed down into large valleys 
and to the plains [23]. These nunataks and unglaciated peripheral 
area could have provided refugia for some species [18] and 
determined several recolonizing routes sculpturing a specific and 
complex genetic structure in the Pyrenees. 

Here we investigate the genetic structure and phylogeographical 
history of the alpine shrub Rhododendron ferrugineum (Ericaceae) in 
the Pyrenees. This species is insect pollinated and reproduces both 
sexually (mixed mating system) and vegetatively through layering. 
It is widely distributed at high elevation, on non-calcareous 
bedrock. The aim of the present study is to understand to what 
extent the phylogeographical history of the species has shaped its 
present genetic patterns. The Asian origin of the genus Rhododen- 
dron [27] implies an east/west migration from Asia to Europe. In 
France, this would lead to the Alps being colonized by R. 
ferrugineum before the Pyrenees, with the Pyrenean populations 
most likely derived from populations from the Alps. After the 
migration to the Pyrenees, the species had to survive successive 
glaciation events. Different hypotheses of glacial refugia and 
postglacial colonization and their consequences on the genetic 
patterns of R. ferrugineum may be proposed: 1 . The species could 
have survived glaciation episodes in one or both of the two main 
refugia at the eastern and south central part of the Pyrenees (HYP 
1; [24]). Thus, the highest genetic diversity should be found in one 
or both of these areas and the genetic pattern would be shaped by 
founder effects i.e. gradual reduction of genomic variability along 
the colonizing routes from main refugia. This would have 
occurred because colonization usually involves only a fraction of 
the genetic diversity present in refugial areas [28]. 2. The species 
could also have survived in nunataks along the mountain chain 
(HYP 2). Hypothetically, high elevation nunatak populations 
would have experienced long periods of isolation, inbreeding and 



genetic drift [18]. This pattern would lead to a patchy distribution 
of genetic clusters with low levels of genetic diversity and high 
population differentiation; 3. A pattern with less differentiated 
clusters than in the previous hypothesis may be observed if the 
species had survived in large lowland refugia in which extensive 
gene flow had occurred. The recolonization of highland areas 
from lowland refugia could have left small, isolated and marginal 
rear-edge populations like those currently observed at low altitude 
along the chain (HYP 3). Indeed, the strong abiotic mountain 
heterogeneity may have allowed remnant populations to match 
suitable ecological conditions at intermediate or low altitudes. 
Moreover, the long life span of the clones [29,30] could have 
permitted their survival even with low seedling recruitment. Such 
rear-edge populations should exhibit more rare alleles and higher 
among-population genetic diversity than populations recently 
established from high altitude source populations [15]. These 
three hypotheses are not exclusive and it is conceivable that the 
genetic pattern observed may result from the different mechanisms 
at work successively or simultaneously within the bioclimatic 
regions of the Pyrenees inducing diverse phylogenetic patterns 
along the chain. 

Materials and Methods 

Species studied 

Sampling of the species of interest was conducted in accordance 
to national and international guidelines (authorization from the 
ONE for the populations BetL, BetH, LecL, PraL, PuiH, from the 
Pare National des Pyrenees for the population TroH, from the 
Gonseil general des Pyrenees-Orientales for population BouH, and 
no permissions were needed for the other sites). 

Rhododendron ferrugineum L. (Ericaceae) is an evergreen shrub with 
a mean height of 70 cm that dominates subalpine landscapes in 
the Alps and Pyrenees between 1500 to 2200 m a.s.l, on north-to- 
west facing slopes, sometimes reaching 90-100% of the vegetation 
cover. It is also found in the Jura Mountains, in the Northern 
Apennines and in the Apuan Alps [31]. In the Pyrenees, 
infrequent smaller and isolated populations can be found lower 
down, between 900—1500 m a.s.l. A single population of some 
scattered individuals remains at low altitude (1500 m a.s.l.) in the 
French Jura Mountains. Rhododendron ferrugineum is a clonal, self- 
compatible [32] mass-flowering species (producing until 3000 
flowers per m^; [33]) pollinated by honeybees and bumblebees 
[32,34]. 

Study sites and sample collection 

We sampled 29 populations across the Pyrenees, the population 
in the Jura Mountains and three populations in the Alps (Figure 1, 
Table 1) ranging from 1070 m to 2200 m a.s.l. The populations 
were divided into three groups according to their elevation: low 
altitude (hereafter L; from 1000 to 1450 m a.s.l), intermediate 
altitude (hereafter I; from 1450 to 1900 m a.s.l.), high altitude 
(hereafter H; from 1900 to 2200 m a.s.l). In the Pyrenees, 
populations were sampled along a longitudinal transect through 
three regions: Languedoc-Roussillon (eastern part), Midi-Pyrenees 
(central part) and Aquitaine (western part). Due to calcareous 
bedrock, populations in the western part of the Pyrenees are rare 
and we only sampled one. An estimation of population size is given 
in Table 1 [<1 ha, 1-10 ha, 10-20 ha, 20-50 ha, 50-100 ha, 
>100 ha). The mean geographical distance between populations 
from Pyrenees vs. Alps, Pyrenees z;^ Jura Mountains and Alps vs. 
Jura Mountains was 900 km (ranging from 515 to 1430 km), 
695 km (from 590 to 830 km) and 335 km (from 165 to 640 km) 
respectively. The distances between populations across the 
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Pyrenees ranged from 3 to 370 km with an average of 126 km. In 
each of the 33 populations (except for the very small JuraL and 
LapL populations), young leaves of twenty individuals were 
sampled (a total of 645 individuals sampled), separated by at least 
5 m in order to avoid sampling clonal individuals, were dried in 
silica gel in the field and then stored at — 20°C. 

DNA extraction and microsatellite markers 

For each individual sampled, 20 mg of leaf tissue was ground 
with two tungsten beads (diameter 3 mm) in a sterile 2 ml 
Eppendorf tube at 30 Hz for 90 s using TissueLyser II (Qiagen). 
Then, the quality and amount of DNA obtained from the DNeasy- 
Plant Mini Kit (Qiagen, Courtaboeuf, France) following the 
manufacturer's instructions was verified using a Nanodrop ND 
1000 spectrophotometer (Peqlab, Erlangen, Germany). 

We used nuclear microsatellites as suitable molecular markers to 
reconstruct regional phylogeographical patterns [18]. Chloroplasts 
(personal observation) or other types of genetic markers such as 
AFLP proved to be little informative in this species [35]. The 27 
polymorphic microsatellite markers (Table SI) were previously 
developed [36,37] using pyrosequencing technologies (454 FLX 
Titanium, Roche Applied Science, Meylan, France). Forward 
primers were labeled with fluorochromes (Eurofins MGW 



Operon, Courtaboeuf, France) and were used in six multiplexed 
PCRs, optimized using Multiplex Manager 1.1 [38]. 

Amplifications were carried out in a 10 |ll mix containing 1 0 ng 
of template DNA, 0.7 x Qiagen Multiplex PGR Master Mix, 
1 0 jiM of each primer and RNase free water. Gycling conditions 
were 15 min at 95°G, 30 x (30 s at 95°G, 90 s at 56°G, 45 s at 
72°G) and 30 min at 60°G. Each set of reactions included a 
negative (water) and a positive (known genotype) control. 
Fragment lengths were read on an ABI 3730 sequencer (Applied 
Biosystems, Gourtaboeuf, France) and scored with Genemapper 
version 4.0 software with the GeneScan - 600 LIZ (Applied 
Biosystems, Gourtaboeuf, France) as internal size standard. 

Data analysis 

Genetic diversity. All 645 individuals were tested for 
clonality with GenAlEx 6.5 software [39]. The presence of null 
alleles at each locus was examined using Micro-Ghecker software 
[40], and tests for departures from Hardy- Weinberg equilibrium 
(HWE) were tested by Fisher's exact tests using Genepop 3.4 [41]. 
We used GenAlEx 6.5 software to calculate for each locus and 
population the number of alleles, the allelic richness (Ar), the 
observed (Hq) and expected heterozygosity (Hg) and the 
inbreeding coeflicient (Fjs). 




Figure 1. Location of the sampled populations, a. Location of the sampled Rhododendron ferrugineum populations in the Alps, Jura Mountains 
and Pyrenees, b. detailed map of the Pyrenean populations (the colors correspond to the five genetic clusters). 
doi:10.1371/journal.pone.0092976.g001 
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Table 1. Geographical information and genetic diversity of Rhododendron ferrugineum populations. 





Moutain 
range 


Sites 


Code 

population 


Location 


Population 
size (ha) 


Altitude 

(m) 


n 


Ar 


Private 
allele 


Ho 


He 


Fis 


Alps 


La Vormaine 


Alp! 


4b Uz U/ N (Jo z4 ly b 


>1 00 


1 560 


20 


3.70 


1 


0.49 


0.49 


0.04 




Col du Lautaret 


AlpH 


46 00 08 N 05 57 18 t 


5-10 


2100 


20 


4.30 


10 


0.59 


0.57 


—0.01 




Parscherkofel 


AutH 


47 12 32 N 1 1 27 18 b 


50-100 


2000 


20 


2.89 


2 


0.41 


0.39 


—0.03 


Jura 




JuraL 


42°57'N 05°44'E 


<1 


1300 


10 


2.67 


0 


0.48 


0.44 


-0.06 


Eastern 
Pyrenees 


Leca 


LecL 


42°28'07"N 02°3r40"E 


<1 


1250 


20 


2.37 


1 


0.34 


0.30 


-0.11 




Prats de Mollo 


PraL 


42°25'21"N 02°24'38"E 


<1 


1450 


20 


2.44 


1 


0.31 


0.32 


0.02 




Lapreste 


LapL 


42°24'34"N 02°22'51"E 


<1 


1230 


15 


2.44 


0 


0.42 


0.34 


-0.21 




Eyne 


EynH 


42°28'27"N 02°07'05"E 


5-10 


1915 


20 


2.78 


0 


0.38 


0.35 


-0.06 




Puigmal 


PuiH 


42°23'26"N 2°05'0"E 


5-10 


2045 


20 


2.74 


0 


0.39 


0.34 


-0.10 




Lac des Bouillouses 


BouH 


42°33'31"N 01°59'42"E 


50-100 


2005 


20 


2.56 


0 


0.37 


0.34 


-0.05 


Central 
Pyrenees 


Camurac 


Cam! 


42°46'30"N 01°55'45"E 


10-20 


1600 


20 


2.92 


0 


0.42 


0.38 


-0.09 




Massat 


Mas! 


42°49'18"N 01°30'75"E 


<1 


1500 


20 


2.59 


0 


0.41 


0.40 


0.0007 




Soulcem 


SouH 


42°39'28"N 01°27'27"E 


50-100 


2035 


20 


2.63 


0 


0.41 


0.39 


-0.01 




Soulcem 


SouL 


42°41'28"N 01°27'07"E 


<1 


1240 


20 


2.74 


0 


0.37 


0.39 


0.08 




Etang de Lhers 


VicH 


42°47'36"N 01°24'16"E 


20-50 


1940 


20 


2.74 


0 


0.38 


0.35 


-0.07 




Etang de Lhers 


VicL 


42°48'23"N 01°22'31"E 


1-5 


1275 


20 


2.44 


0 


0.36 


0.32 


-0.09 




Espot 


Spol 


42°33'47"N 01°05'36"E 


<1 


1515 


20 


4.18 


0 


0.35 


0.43 


0.10 




Bethmale 


BetL 


42°51'42"N 01°04'02"E 


1-5 


1385 


20 


2.67 


0 


0.40 


0.37 


—0.06 




Bethmale 


BetH 


42°50'24"N 01°03'42"E 


50-100 


1920 


20 


2.92 


0 


0.40 


0.37 


—0.03 




Val d'Aran 


AraH 


42°39'45"N 00°58'56"E 


<1 


2075 


20 


2.26 


0 


0.29 


0.29 


0.05 




Le Mourtis 


MouL 


42°54'32"N 00°46'31"E 


5-10 


1440 


20 


2.59 


0 


0.33 


0.33 


0.04 




Peyragudes 


GudI 


42°47'42"N 00°27'14"E 


5-10 


1650 


20 


2.67 


0 


0.39 


0.37 


-0.03 




Hourquette 
d'ancizan 


Houl 


42°53'57"N 00°18'13"E 


5-10 


1565 


20 


2.74 


0 


0.42 


0.40 


-0.04 




Neouvielle 


NeoH 


42°48'N 00°09'E 


50-100 


2200 


20 


2.56 


0 


0.47 


0.44 


—0.03 




Le Peyras 


PeyL 


42°59'15"N 00°10'46"E 


5-10 


1220 


20 


2.33 


0 


0.36 


0.36 


0.03 




Col du Tourmalet 


TouH 


42°53'57"N 00°07'54"E 


5-10 


1950 


20 


2.81 


1 


0.42 


0.39 


-0.04 




Cirque de 
Troumouse 


TroH 


42°43'41"N 00°05'40"E 


5-10 


2080 


20 


2.44 


0 


0.39 


0.38 


0.02 




Chiroulet 


ChiL 


42°57'38"N 00°04'17"E 


1-5 


1240 


20 


2.70 


0 


0.38 


0.37 


0.02 




Luz-Ardiden 


LuzH 


42°52'45"N 00°03'34"W 


10-20 


1905 


20 


2.63 


0 


0.37 


0.37 


0.02 




Huesca 


Espl 


42°45'08"N 00°14'18"W 


>100 


1640 


20 


2.27 


0 


0.20 


0.30 


0.32 




Lac du Tech 


TecL 


42°54'40"N 00°15'32"W 


<1 


1240 


20 


2.41 


0 


0.39 


0.34 


-0.11 




Col du Somport 


SomI 


42°47'4V'N 00°31'38"W 


<1 


1645 


20 


2.26 


0 


0.33 


0.37 


0.16 


Western 
Pyrenees 


Col de la Pierre 
Saint Martin 


Pie! 


42°58'07"N 00°45'52"W 


<1 


1830 


20 


2.00 


0 


0.28 


0.29 


0.09 



Altitudinal classes are given in the code population (L: low, I: intermediate, H: high altitude). 

n = sample size, Ar = allelic richness. Ho = observed heterozigosity. He = expected heterozygosity, Fis 

values are in bold, P<0.01). 

doi:l 0.1 371/journal.pone.0092976.t001 



within population coefficient of inbreeding (non-significant 



Population structure and diflerentiation. The pairwise 
FsT (level of genetic differentiation; [42]) among populations were 
calculated with Genodive v2.0 [43]. Then, the statistical signifi- 
cance of the pairwise Fgx values was assessed by 1000 random 
permutations of individuals across populations with Genetix 4.05 
[44]. 

To analyze the genetic structure among populations, we used 
Structure software (ver. 2.3.3; [45]) in order to determine the most 



likely number of homogenous clusters in the sample (K). The 
admixture model was applied with 50 runs for each K value from 
1 to "number of populations +1", each run comprising a burn-in 
period of 20 000 iterations followed by 100 000 iterations. The 
optimal value of was evaluated by considering the highest mean 
likelihood value of A', i.e., L (A), as well as the AA^ method [46]. 
These results were computed and visualized in the online interface 
Structure Harvester [47]. 
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To evaluate the distribution of genetic variance within and 
among populations, and among groups when populations were 
nested into genetic clusters, altitude (low, intermediate and high 
altitude) or population size (according to the six population size 
groups in Table 1) we performed analyses of molecular variance 
(AMOVA; [48]) using Genodive v2.0. The significance of the 
results was tested using 999 random permutations of the data. 

Mantel tests [49] were performed ("ncf ' package implemented 
in R, R Development Core Team, 2008) to test for isolation-by- 
distance (IBD; divergence due to drift and mutation proportionally 
increasing with geographical distance), investigating the correla- 
tion between matrices of genetic (i^sx) and geographical (km) 
distances. Mantel test were also used to determine wheter 
elevational separation was related to genetic differentiation. 
Calculation of the gene flow on the basis of rare alleles was 
conducted [50] and the number of migrants (N^) was corrected for 
the number of samples. 

Results 

A total of 153 alleles at 27 microsatellite loci were detected 
across 645 individuals at all sites (Table S2). No clones were 
identified among the individuals sampled as revealed by the test of 
clonality. 

Significant departures from HWE across loci were detected in 
12 out of 33 (or 36%) sites. The RF87P2 marker was the only locus 
that showed no departures from HWE. Most of the loci showed 
significant departures from HWE in two to seven populations, 
except RF14P3, RF47P1 (11 populations each), RF56P1 (13 
populations) and RF96P2 (31 populations). These instances likely 
reflect occasional departures from random mating rather than the 
presence of null alleles. Indeed, the probability of null alleles was 
low for all loci, according to Oosterhout's method [40], mean 
estimation of null alleles was —0.031 (ranging from —0.3675 to 
0.1352). No differences in the results were observed when 
including the microsatellite marker RF96P2 or not, so we decided 
to keep it in the analyses. 

Genetic diversity and allelic richness among populations 

The number of alleles ranged from 2 to 15 (mean 6.25) alleles 
per locus. A total of 16 private alleles (restricted to a single 
population; Table 1) was detected (10.45%), mostly in AlpH (10), 
AutH (2), AlpI (1) and also in some Pyrenean populations (TouH, 
PraL, LapL had one private allele each). Allelic richness (Ar) 
varied across sites (2-4.30). The highest Ar were observed in the 
Alps (AlpI = 3.70 and AlpH = 4.30) and in one of the southernmost 
populations of the central Pyrenees (Spol = 4.18), while the lowest 
Ar were observed at the western end of the Pyrenees (Piel = 2.00, 
SomI = 2.26, Espl = 2.27) and in the central Pyrenees 
(AraH = 2.26). Ar, Hq and He were 4.6, 0.37 and 0.35 
respectively in the Pyrenees, and 5.2, 0.5 and 0.45 in the Alps. 
Expected (He) and observed heterozygosity (Hq) ranged from 0.20 
to 0.59 and from 0.29 to 0.56 respectively for all sites. He was high 
in the Alps and Jura Mountains (ranging from 0.39 to 0.57), and in 
some Pyrenean populations (Spol and NeoH, 0.43 and 0.44 
respectively) and the lowest in the easternmost (LecL, PraL, AraH) 
and westernmost (Espl, Piel) Pyrenean populations. These latter 
populations were small and isolated (< 1 ha) except the large and 
central population Espl (>100 ha). Significant heterozygosity 
excesses (Fjs coefficients ranging from —0.07 to —0.21) were 
detected in three eastern (LecL, LapL, PuiH) and four central 
populations of the Pyrenees (CamI, VicH, VicL, TecL). In 
contrast, significant heterozygosity deficiency was found in four 
populations with particularly high values of Fjs coefficients in 



SomI (0.16) and Espl (0.31). Because the frequencies of null alleles 
did not appear to be significant, the highest Fjs values could be 
due to inbreeding. 

Genetic differentiation and structuration among 
populations 

The Est values (Table S3) were significantly different from zero 
for all population pairs, except for three pairwise comparisons: 
BetH/BetL and SouH/SouL that were geographically really close 
(less than 3 km) and Caml/EynH (35 km apart). The mean Est 
value was relatively low (0.267) and most of the highest values 
(>0.33) were found in the Pyrenees between eastern and western 
populations, reaching 0.45 for LecL/ AraH. The population from 
the Jura Mountains showed a highly variable level of Egx (from 
0.25 to 0.42). 

To infer the relationship between populations we used the 
clustering algorithm implemented in the program Structure. 
Comparing all 33 populations (K: 1 to 34), AK as a function of 
K reached a peak at K = 2 and then reached a second peak at 
K = 6 (Figure 2a). At K = 2, populations from the Alps, Jura and 
the eastern Pyrenees were assigned to the same cluster while 
populations from western and central Pyrenees were clearly 
grouped apart. Individuals from AlpI showed a similar proportion 
of membership to both genetic clusters. At K = 6, populations from 
the Alps and Jura were assigned to one cluster (CLl) and those of 
the Pyrenees were grouped in 5 clusters (CL2 to CL6) along a 
longitudinal gradient. Except between CL5 and CL4 or CL6, 
there was low admixture between the other Pyrenean clusters. 

The analysis of molecular variance (AMOVA; Table 2) based 
on Fgx showed that as expected for polymorphic loci such as 
microsateUites, the within-population variance component was 
rather high (>66%, P = 0.001). The among variation component 
(Est) ranged from 9% to 31% according to the grouping and was 
statistically significant in all combinations tested (P = 0.001). At a 
superior hierarchical level, the genetic clusters could explain a 
significant part of the total genetic variance (24.4%, P = 0.001) 
whereas it was not significant when grouping populations into 
altitudinal groups (P = 0.956) or population size (P = 0.222). 

Allelic richness (Ar) was highest in CLl (Alps and Jura, 5.2) and 
lowest in CL6 (western Pyrenees, 2.9) and ranged from 3.2 to 3.8 
in the other clusters (Table 3). Fsx values were all significantly 
different among clusters (Table 4). The highest Est values were 
found between CL2 (eastern Pyrenees) and CL5 or CL6 (western 
Pyrenees, 0.29 and 0.31 respectively). The Est values between 
CLl and the Pyrenean clusters were quite low (ranging from 0.16 
to 0.20). 

There was a positive relationship (Mantel test) between genetic 
distance and geographical distance matrices across populations 
located in the Pyrenees (R = 0.784, P = 0.001). The increase of 
isolation with geographical distance was consistent with the 
previous results (see Structure analyses). When including popula- 
tions from the Alps and Jura, the Mantel test was still positive 
(R = 0.26) but not significant (P = 0.383), confirming the lack of 
strong differentiation among the populations from the Alps and 
Pyrenees. At the Pyrenean scale, the correlation among elevational 
distances and EST was slightly positive but not significant 
(R = 0.092, P = 0.06). 

Gene flow among populations 

Gene ffow among all studied populations (Alps, Pyrenees and 
Jura) was low and reached a mean value of Nm = 0.72 migrants 
per generation. The values reached 2.88 for all Pyrenean 
populations but 2.10 in the eastern lineage (CL2) and 5.9 in the 
western lineage (CL3 to CL6). 
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Figure 2. Genetic structure of the 33 populations, a. Plot of Delta K according to K. b. Structure clustering results obtained at K = 2 and K = 6. 

Each individual is represented by a thin bar corresponding to the sum of assignment probabilities to the K cluster. Black bars separate populations. 
doi:10.1371/journal.pone.0092976.g002 



Discussion 

Genetic diversity in R. ferrugineum populations 

Microsatellite analysis revealed a relatively high level of genetic 
variability in R. ferrugineum populations in the Pyrenees (Ar: 2-4.18, 
He: 0.29-0.44) but lower than in the Alps and the Apennine 
populations (Ar: 2.71-5.57, He: 0.37-0.68; [31]) or in other long- 
lived species found in the Pyrenees {Pinus uncinata: He: 0.83-0.99, 



Table 2. Partitioning of the observed microsatellite variation 
(AMOVA) based on fsj of Rhododendron ferrugineum 
populations. 





Source of 
Variation 


SSD 


d.f. 


MS 


Var- 
comp 


% 

Variation 


P- 

value 


Within Population 


4355.120 


548 


7.947 


7.947 


0.662 


0.001 


Among Population 


732.618 


24 


30.526 


1.135 


0.095 


0.001 


Among Genetic 
cluster 


1461.242 


4 


365.310 


2.929 


0.244 


0.001 


Within Population 


4354.038 


548 


7.945 


7.945 


0.690 


0.001 


Among Population 


2103.864 


26 


80.918 


3.659 


0.310 


0.001 


Among Altitude 


93.988 


2 


46.994 


-0.167 


0.00 


0.956 


Within Population 


4366.874 


548 


7.894 


7.894 


0.685 


0.001 


Among Population 


2079.120 


23 


83.837 


3.372 


0.293 


0.001 


Among Population 
size 


106.877 


5 


56.922 


0.258 


0.022 


0.222 
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(25); Quercus petraea: He: 0.76-0.84, [51]). Although clonality and 
selfmg are known to decrease genetic variability [52], the long life 
span of individuals [29,30] and the strong life-time inbreeding 
depression eliminating most selfed individuals (personal data) may 
contribute to maintaining relatively high genetic diversity even in 
small R. ferrugineum populations. Accordingly, most of the 
populations studied showed a Fjs close to zero and a great 
proportion of the genetic variation was found within populations 
as observed in other R. ferrugineum populations [31] and other 
alpine species [19], [53], [18]. However, selfmg has probably 
contributed to individual recruitment in populations harboring 
higher heterozygote deficiency (SomI, Piel, Espl, Spol). Such high 
heterozygote deficiency has also been reported in the Apennines 
and the Alps [31]. Surprisingly, significant heterozygote excess 

Table 3. Genetic diversity of Rhododendron ferrugineum per 
genetic cluster. 





Ar 


Ho 


He 


F,s 


Cluster 1 


5.630 


0.498 


0.584 


0.148 


Cluster 2 


3.630 


0.375 


0.364 


-0.030 


Cluster 3 


3.185 


0.386 


0.406 


0.048 


Cluster 4 


3.778 


0.360 


0.390 


0.078 


Cluster 5 


3.481 


0.402 


0.429 


0.062 


Cluster 6 


2.926 


0.315 


0.384 


0.179 



Ar = allelic richness, Hq = observed heterozigosity, He = expected 
heterozygosity, F|s = within population coefficient of inbreeding. 
doi:10.1371/journal.pone.0092976.t003 
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Table 4. fsj values among genetic clusters (p-values for all pairs<0.001). 







Cluster 1 


Cluster 2 


Cluster 3 


Cluster 4 


Cluster 5 


Cluster 6 


Cluster 1 




0.203 


0.171 


0.177 


0.165 


0.166 


Cluster 2 






0.237 


0.275 


0.297 


0.311 


Cluster 3 








0.101 


0.153 


0.230 


Cluster 4 










0.108 


0.235 


Cluster 5 












0.129 


Cluster 6 
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occurred mostly in eastern (LapL, LecL, PuiH), and central 
populations (CamI, VicH, VicL) and in one of the westernmost 
populations (TecL). 

Alpine origin of the two Pyrenean R. ferrugineum lineages 

The Rhododendron genus is known to originate from Asia [27] and 
colonized Europe from east to west. Our hypothesis was that 
migration from the Alps to the Pyrenees occurred with a loss of 
genetic diversity along the colonization routes. In our study, we 
identified two lineages in the Alps: the first, found only in the 
northern population (Alp I), was related to western Pyrenean 
populations. This lineage would possibly be found in more 
populations with more extensive sampling in the Alps. The second 
lineage, related to the eastern and central Pyrenees, was found in 
the same population AlpI but also in south and central-eastern 
populations from the Alps and in the Jura Mountains. At K = 6, 
populations from the Jura Mountains and the Alps are found in 
the same cluster, that may indicate the occurrence of recent gene 
flow between the Alps and the Jura Mountains, possibly in low 
land refugia during glaciations. It would be interesting to complete 
our data by a large sampling of populations from the Alps and the 
Jura Mountains in order to detail the European phylogenetic 
pattern of the species. Compared to the Alpine populations, 
Pyrenean populations were characterized by a substantial loss of 
genetic diversity (Ar, Hg, Hq) suggesting bottlenecks during the 
migration events [54], [16]. Cold, dry steppe /tundra occupied 
most of the land between the Alps and the Pyrenees during the last 
maximum glaciation [55] and could have allowed migration of 
Alpine species such as R. ferrugineum either from the southwestern 
Alps [56], [57] or from northern Alpine populations through the 
French Massif Central which is considered to be a key area in the 
colonization of the Pyrenees by Alpine species [58] . 

Genetic pattern and phylogeographical history of 
/?. ferrugineum in the Pyrenees 

We found that genetic diversity of R. ferrugineum was not 
correlated to altitude or population size; instead there was a sharp 
longitudinal split between the two genetic lineages. Other studies 
in alpine ecosystems similarly found no changes in genetic diversity 
depending on population size or along altitudinal gradients [59]. 
Along environmental gradients, several parameters such as overlap 
of generations or recruitment frequency may influence intrapop- 
ulation genetic variation and gene flow among populations may be 
sufficient to maintain relatively high genetic diversity even in the 
smallest populations. 

The longitudinal genetic pattern seemed to result from the 
migration of two Alpine lineages and from the establishment of 
populations in refugia along the Pyrenean chain due to glaciation 
events. The strong differentiation between the eastern and western 



Pyrenean lineages and the low level of gene flow observed between 
the two suggested that very few admixtures occurred since the two 
lineages had become established. This may result from topo- 
graphical barriers (high summits, glaciated areas), seed dispersion 
by gravity and short pollinator travels hampering gene flow among 
populations and genetic homogenization. An east-west genetic 
separation had also been observed in other Pyrenean species (i.e. 
Papaver alpinum, [57], Tr folium alpinum, [19]). The populations of 
the eastern lineage were little differentiated even though they were 
distant from one another (35 km between the two undifferentiated 
populations CamI and EynH), were of very different sizes, grew in 
contrasting environments, and experienced a low level of among- 
population current gene flow. Very small and isolated populations 
(PraL, LecL, LapL) showed a relatively high level of intrapopu- 
lation heterozygosity (H^, Hq) and exhibited private alleles (PraL, 
LapL) despite relatively lower allelic richness likely due to the loss 
of rare alleles by genetic drift, common in such small populations 
[60]. These populations survive in rather unusual ecological 
conditions for this species, i.e. on north facing slopes in deep, 
shaded valleys and mainly beneath beech forest {Fagus sylvaticd) or 
sometimes in the vicinity of Mediterranean Quercus ilex forests. 
These data suggested that they were more likely to be rear-edge 
remnant populations [12] left there after the recolonization of 
highland areas (HYP 3) rather than the result of a secondary 
colonization of lowland areas by high altitude individuals. 
Together, these findings and the absence of clear founder effects 
suggested that the eastern populations within CL2 have experi- 
enced recent and extensive gene flow. Such intense gene 
exchanges could have occurred in large lowland refugia during 
the cold periods of the Pleistocene and possibly later in genetically 
connected populations growing over high altitude plateaus, which 
are typical of eastern Pyrenees landscapes. These large high 
altitude populations could have then become fragmented possibly 
during the last 30 ky when climate became warmer and drier [26]. 

The four well separated and differentiated genetic clusters 
nested in the western Pyrenean lineage appeared as a sub-structure 
of the east/west configuration and presumably diverged more 
recently. Overall, we did not observe a gradient of genetic 
variability loss among these clusters as expected from a range- 
expansion wave [11]. Rather, we observed a consistent genetic 
heterozygosity and allelic richness and a significant isolation by 
distance (Mantel test) among them. Wright's isolation by distance 
model suggests that gene flow is reduced as geographical distance 
increases. This may explain the very low level of admixture 
between CL3/CL6. The patchy distribution of often homogenous 
and highly distinct clusters enclosing a similar level of within- and 
between-population genomic variability is compatible with the 
hypothesis of the survival of the species in numerous Pyrenean 
refugia as proposed for other Pyrenean species [19], [18]. 
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Moreover, the level of differentiation, the sharp genetic boundaries 
between most clusters and the number of migrants per generation 
suggested a long isolation of populations and the action of genetic 
drift in several high altitude refugia (HYP 2) rather than the 
survival in large and genetically connected populations in low 
elevation unglaciated areas. The higher genetic diversity in Spol 
likely accounts for such refugia as well as TouH (having one 
private allele) where genetic variability was greater than in other 
populations. However, the relatively higher level of admixture 
between CL5 and CL4 or CL6 suggested that they had recently 
exchanged genes. Because the admixture concerned mainly 
isolated and sometimes very distant (AraH and CL5 populations) 
populations at high and intermediate altitude, contemporary gene 
exchanges among them through pollen or seed dispersal was very 
unlikely. We suspected that, at least in this part of the Pyrenees, 
longitudinal gene exchanges could have occurred possibly in 
lowland refugia across Pyrenean piedmonts before the post-glacial 
recolonization of high altitude area through valleys (HYP 1). Such 
a recolonizing pattern has been observed in sessile oak [51] and 
suggested for Aster pyrenaeus [20]. Thus, the phylogeographical 
pattern in central and western Pyrenees could result both from 
population isolation and differentiation in nunataks and gene 
exchange in some lowland refugia. Among clusters, CL6 showed 
the lowest genetic variability (Espl, SomI and Piel) and genetic 
diversity gradually decreasing in this cluster along the sequence 
TecL, SomI and Piel. Thus, the genetic pattern of CL6 could be 
chiefly related to founder events and the loss of genetic diversity 
and heterozygosity [1 1] during westward expansion of the species. 

Implications for conservation 

The longitudinal genetic structuration found in R. ferrugineum 
and particularly the east/ west separation between the two main 
genetic lineages has already been observed in other studies [57], 
[19]. This suggests that many other alpine species may have been 
subject to similar events as R. ferrugineum and, nowadays, could 
present a comparable genetic pattern. We state that these data 
now should be usefully considered in management projects for 
Pyrenean species particularly those having an extended east/ west 
distribution. According to the European Economic Community 
laws (92/43/EEC habitat directives, 1992) R. ferrugineum heath- 
lands must be locally preserved due to their high ecological and 
patrimonial value. We think that, henceforth, in aim to preserve 
genetic diversity of the species, the genetic pattern here revealed 
should be taken into consideration to choose sites that have to be 
protected. The eastern lowland isolated populations (LecL, PraL, 
LapL) grow in very specific habitat and present rare alleles. 
According to the "rear edge population hypothesis", they harbor 
all together the bulk of species' genetic diversity [12]. Scientific 
programs have to be developed to follow the population dynamic 
and the genetic diversity evolution of these very original 
populations and take an inventory of other currently unknown 
rear-edge populations. Moreover, studying population dynamics at 
the edge of a species' range is crucial to understand the response of 
species to environmental changes and their susceptibility to 
extinction [31]. Also the populations TouH (presenting a rare 
allele) and Spol (having a high allelic richness) should be 
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